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ABSTRACT 

Recent developments of new medical treatment techniques put challenging demands on ultrasound imaging 
systems in terms of both image quality and raw data size. Traditional sampling methods result in very large 
amounts of data, thus, increasing demands on processing hardware and limiting the flexibility in the post- 
processing stages. 

In this paper, we apply Compressed Sensing (CS) techniques to analog ultrasound signals, following the recently 
developed Xampling framework. The result is a system with significantly reduced sampling rates which, in turn, 
means significantly reduced data size while maintaining the quality of the resulting images. 
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1. INTRODUCTION 

Modern ultrasound systems utilize an array of transducer elements in a process known as beamforming 00 An 
imaging cycle begins when modulated acoustic pulses are transmitted from some or all array elements. Specific 
delays are applied to the transmitting elements, such that the interfering waves form a narrow beam, along which 
most energy propagates. As the energy propagates, the beam gets narrower, until reaching the focal point, after 
which it expands. Determining the position of the focal point is achieved by applying appropriate delays to the 
transmitting elements. 

As the focused energy propagates along the tissue, echoes are scattered and reflected by density and propagation 
velocity perturbationsP The array elements detect the reflected energy. A second beamforming process is 
then performed, aimed at localizing the reflecting elements, while at the same time improving signal to noise 
ratio (SNR). This beamforming process is performed by summing delayed samples of the received data. For that 
purpose, high rate A/D conversion must be first carried out at each of the receiving channels. Actually, optimizing 
the image resolution requires that each receiving channel will be sampled at 3-5 times the center frequency of 
the modulated pulseP More explicitly, let us consider a B-Mode scan, in which the central frequency of the 
transducer may vary in between 2 — 15mhz, depending on the useP Advanced composite materials, often used 
in the transducer, can attain a relative bandwidth in excess of 100%. Therefore, if we assume a nominal center 
frequency of 5m Hz, we end up with a baseband bandwidth of approximately 10m Hz- Confined to the classic 
Nyquist-Shannon sampling theoremp where the only prior on the signal is that it is bandlimited, this implies 
that standard ultrasound devices must sample the analog signal received in each active element at a typical rate 
of at least 20mhz- 

Recent works^ED show that, by exploiting other priors regarding the signal structure, it is possible to design more 
efficient sampling schemes, which actually break through the Nyquist barrier. Decreasing the sampling rate is of 
vast interest, as it may obviously be manifested in reduction of machinery size and power consumption. 
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In their recent work, Tur, Eldar and Friedman^! first proposed to implement these ideas in ultrasound imaging. 
They suggested to regard the signal received in each of the array transducer elements as having finite rate of 
innovation (FRI) P More specifically, they assume each such signal comprises at most L replicas of a known-shape 
pulse, all received within the time interval [0, r). Delay and gain parameters are associated with each replica, 
such that the entire problem may be characterized by 2L degrees of freedom. They then developed a new sub- 
Nyquist sampling technique that exploits this structure, in order to reduce the sampling rate way beyond that 
which is used in standard ultrasound devices, dictated by the classic Nyquist-Shannon sampling theorem. 

Tur, Eldar and Friedman's work follows the spirit of analog compressed sensing, also referred to as Xampling. 6 
The latter is an emerging framework, which combines classic methods from sampling theory together with recent 
developments in compressed sensing, aimed at sampling analog signals far below the Nyquist rate. Throughout 
this paper, we will use the term Xampling whenever referring to the sub-Nyquist sampling schemes purposed in 
Ref. [ToJ and in later work by Gedalyahu, Tur and EldarP^ It should be noted, though, that both these schemes 
are special cases of Xampling, in which the signal's FRI property is exploited in order to achieve the goal of 
sub-Nyquist sampling. 

The Xampling scheme suggested by Tur, Eldar and Friedman^ comprises the following steps: first, the received 
signal is filtered using the compactly supported Sum of Sines kernel. The filtered signal is then sampled at 
nearly the rate of innovation, namely the number of unknown parameters 2L, which is much smaller than the 
Nyquist rate of the pulse. The extracted samples are used for computing a finite set of Fourier coefficients, 
which correspond to the r-periodic extension of the received signal. Having obtained the finite set of Fourier 
coefficients, spectral analysis techniques, similar to these presented by Vetterli et al.JSElare utilized, in order to 
estimate the set of unknown delays and amplitudes characterizing the reflected pulses. 

Whereas Tur, Eldar and Friedman^ adopt a filtering and sampling approach as a preliminary step for obtaining 
the required set of Fourier coefficients, Gedalyahu, Tur and Eldar'i^ later work suggests a second approach for 
obtaining the same coefficients: the single filtering and sampling channel is replaced by a bank of modulators 
and integrators. Both approaches are aimed at minimal rate sampling of a single received channel. 

Referring to an array of transducer elements utilized in the ultrasound imaging device, by processing the signal 
received in each channel separately (using either approach), one may obtain a corresponding set of delays and 
amplitudes from low-rate samples. The sets obtained from all receiving elements may then be combined (via 
some geometric interpretation), in order to estimate the two-dimensional coordinates of the reflecting elements. 
Note, however, that such schemes cannot achieve the SNR improvement which is an integral part of standard 
beamforming techniques; this is because the correlation between signals received in different channels is not 
exploited throughout the process of extracting the parameters (pulse amplitudes and delays) from each signal. 
More explicitly, the Xampling schemes proposed in Rcfs. [10| and \TT\ both aim at accurately detecting strong, 
localized pulses, related with macroscopic perturbations. However, the actual signals received by the array 
elements also contain components which, in the context of our work, may be regarded as noise. Actually, part 
of these noisy components arises from constructive and destructive interference of acoustic waves reflected from 
dense, subwavelength scatterers in the tissue (these are typically manifested as granular texture in the ultrasound 
image, called speckle, after a similar effect in laser opticfl). Apparently, these noisy components induce erroneous 
results when Xampling the received signals. If we wish to obtain meaningful results by either Xampling scheme, 
while maintaining a rather low oversampling factor, the overall SNR improvement is indeed a crucial step. 



This paper is aimed at generalizing the schemes of Refs. II and 10 to multiple antenna arrays. Our goal is 
to obtain a two-dimensional, focused ultrasound image, corresponding to strong perturbations in the scanned 
plane, while reducing the sampling rate in each active element, by a factor of 10-15 times relative to the rate 
used in standard ultrasound devices. Furthermore, we aim at achieving this goal in the presence of noise in 
the received signals. In such case, straightforward implementation of either schemes of Refs. [10] or [TT] on each 
channel independently, would require hard thresholding (which in turn cancels/attenuates desired pulses) and/or 
increasing the oversampling factor, such that the final sampling rate grows towards the Nyquist rate. 

Our Xampling scheme's most expensive computational component regards the extraction of the pulses' delays 
and amplitudes from the set of low rate samples. Referring to standard ultrasound devices, this component 
substitutes both the Hilbcrt transform, which is applied throughout the process of envelope detection^ and 



the expensive beamforming computations, namely: summing samples obtained at Nyquist rate from all active 



elements. In addition, referring to straightforward implementation of either systems of Refs. 10 or 11 on each 
channel separately, our scheme extracts the pulses' parameters once per image line, rather than once per active 
element participating in the image line generation (which equivalently means tens of times per image line, 
depending on the number of active elements used for beamforming) . 

The Xampling system we propose may be summarized as follows: Let us assume that we could somehow generate 
the beamformed signal, corresponding to a single image line, in the analog domain. Assuming that such a signal 
maintains the FRI property of the received signals from which it was constructed, we may Xample it using the 
scheme suggested by either Ref.[lO]or[TTJ yielding the delays and amplitudes of pulses along the final beamformed 
image line. By mathematically formulating the process of Xampling the beamformed signal using the scheme 
suggested in Ref. |11[ and then applying several algebraic manipulations to this formulation, we end up with 
a new set of generalized modulation kernels. The latter may now be applied directly to the analog signals 
which are received in each of the active elements, thus bypassing the impractical step of actually generating the 
beamformed signal in the analog domain. 

Due to space limitation, the following paper is not aimed at presenting a rigorous review of our results. Instead, 
we present an outline of our approach, and preliminary results obtained using actual ultrasound data. The 
paper is organized as follows: Section [2] briefly outlines principles of standard ultrasound imaging, namely the 
process of beamforming in both Polar and Linear scan methods, applying dynamic receive (Rx) focus. Section [3] 
reviews the one-dimensional Xampling scheme suggested by Refs. [10] and [TT] In section [4] we present our system, 
which combines the concepts of beamforming with those of one-dimensional Xampling. Section [5] provides results 
obtained by applying our Xampling scheme upon actual ultrasound data. In this section we also analyze the 
reduction in the necessary amount of samples, and the way this affects the overall computational cost. Finally, 
conclusions are drawn in Section [U 

2. DYNAMIC FOCUSING IN POLAR AND LINEAR SCAN 

This section is aimed at outlining the method by which an ultrasound image is generated using a linear array 
of transducer elements. Our discussion refers mainly to B-mode scan, in which the array simultaneously scans a 
plane through the body, resulting in a two-dimensional image, which may be viewed on screen. We are specifically 
interested in the process carried out by the electronic beamformer, in which multiple signals, received from a set 
of receivers, are focused into a single trace, known as the beamformed signal. The latter is envelope detected, 
forming a single image line. The analysis reviewed throughout this section is based mainly on Ref. [3j 

An ultrasound image consists of roughly 100 lines. Each line is generated throughout a single transmit-receive 
cycle. In such a cycle, a set of active elements first transmits acoustic pulses, modulated to central frequency of 2 
to 15m Hz, depending on the use. An appropriate delay is applied to the pulse transmitted from each transducer 
element, aimed at obtaining constructive interference at a specific point in the plane, referred to as the focal 
point. The interfering acoustic waves form an acoustic pulse, which propagates along a narrow beam (containing 
most of the transmitted energy). The beam gets narrower, until reaching the focal point, after which it expands. 

The velocity at which the pulse propagates will be denoted by c, and varies between 1446 m / sec (fat) to 1566 m / sec 
(spleen). An average value of 1540 m / sec is assumed by scanners for processing purposes. As the transmitted pulse 
propagates inside the tissue, it encounters density and propagation- velocity perturbations. These cause scattered 
and reflected echoes, which are detected by the array elements. Applying the acoustic reciprocity theorem,^ a 
second beamforming process is now carried out, in which the received signals are combined into a single trace, 
which in a sense, visualizes structures in the tissue, along the transmitted beam. Modern ultrasound devices 
typically handle 64 to 192 transducer elements in the beamforming process. We emphasize, that combining 
the received signals is performed in the digital domain, implying that modern ultrasound devices must first 
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sample the signal received in each of the active elements at the Nyquist rate (typically 20m Hz)- Refs. 10 and 
introduced an approach, which allowed reconstruction of the analog signal detected by an individual transducer 
element from a very low number of samples. In contrast, our generalized scheme is aimed at obtaining the low 
rate samples from all active elements, in a manner which will enable to directly reconstruct the beamformed 
signal. This is further discussed in Section |4j 



We would now like to formulate the manner in which the beamformer combines the signals received in all active 
elements into a single trace, known as the beamformed signal, and better understand the significance of the 
latter. This formulation applies to conventional ultrasound imaging, and will be necessary when we translate 
our theoretical scheme of Xampling the beamformed signal, into an applicable scheme where the signals received 
in the active transducer elements are sampled directly. 

Referring to Figure [l] we examine the two-dimensional plane XZ in which 2M + 1 elements are aligned along 
the x axis (the center of the array coincides with the origin). We analyze one cycle, in which a single image line 
is constructed. The image line corresponds to a beam located within the XZ plane, emerging from the array 
center. We denote by a, the angle between the beam and z axis (normal to the array). 

The cycle begins, when each active element transmits a single modulated pulse, such that the interference pattern 
may be observed as a concentrated pulse of energy, propagating along the beam. We regard the pulse as if it 
was transmitted from the array center, at a known time, which we shall define as t — 0. Knowing the speed in 
which the pulse propagates (denoted by c, and assumed 1540 m / sec ), we may now estimate the distance which 
it traveled along the beam by the time instance t n , denoted by r(t n ), and thereby its two-dimensional position, 
P n : 

Pn = [ cosa sina ] r (t n ) = [ cosa sina ] ct n . (1) 

Assuming that an echo was reflected due to some perturbation located at p n , we may easily estimate the time 
in which it will arrive back at the origin, To (t n ): 

T (t„) =t n + -||p n || = 2t n . (2) 



Since the distance from p n to each of the active elements varies, each element will detect the reflected pulse at 
a different time instance. Namely, the time at which the pulse arrives at the mth receiver, positioned at x m , is 
given by: 

1 1 / ~ 2 2 

T m (t n ) = t„ + - Pn - X m = t„ + -A/ (ct n Sma - S m ) + (ct n COSa) , (3) 
c c v 

where S m denotes the x coordinate of the receiving element. 

Let us denote by (p m (t), the analog signal detected by the active element indexed m. Beamforming is now achieved 
by shifting each of the received signals ip m (t), in order to compensate for the time difference T m (t n ) — To(t n ), and 
then summing the shifted versions. The acoustic reciprocity theorem^ implies that when we sum the shifted 
signals, constructive interference will occur at t = To(t n ) = 2i n , providing that an echo was indeed reflected from 
p n at time t n . Denoting by $(£;cn) the sum of the delayed signals, we are thus interested in the value which 
<!?(£; a) obtains at To(t n ) = 2t n . 

Shifting the signal <p m (t) so that the difference T m (t n ) — To(t n ) is compensated, is obtained by applying the 
(possibly negative) delay 

m [t n ] a) = T (t n ) - T m (t n ) = t n - \Jt 2 n + {S m /cf - 2t n (Sm/c) sina (4) 

to the signal received in the mth element. 
Summarizing the above, we have: 

M 

$(2t n ;a) = $(£;a)| t =2t„ = <P™ {t - m {t n ;a))\ t =2t n 

m=— M 



m=-M 



= X! Vrn[t„ + ^Jtl + (S m /c) 2 - 2t n (S m /c) sina , (5) 



m - -M 




Figure 1: 2M + 1 elements aligned along the x axis. The ray along which the pulse propagates forms an angle 
a with z axis. We analyze a reflection emerging from the point p n , which is associated with the nth focal zone. 



which corresponds to the intensity of a reflection originating at time t n , from the coordinate p n . For purposes 
of convenience, we will finally substitute 2t n — > t, obtaining an expression for the beamformed signal: 

M ( 1 \ 

$(*; a) = frn ( 2 (* + Vt 2 + 4 (6 m /c) ((<J m /c) - feina)) J . (6) 

m=-M ^ ' 

Referring to ([7J), our last substitution implies that <£>(i;a) represents the intensity of a reflection originating at 
time t/2 from a point distanced ct/2 from the origin, along the transmitted beam. The assumption is that the 
transmitted pulse indeed intersected this point at t/2, and was possibly scattered. a) is obviously the result 
of varying the receive focal point along time, and its construction is therefore referred to as "Dynamic Focusing" . 

We note that, instead of obtaining the beamformed signal using dynamic focusing process, beamformers some- 
times divide the image line into N segments, called focal zones, such that an entire segment of $(t;a), corre- 
sponding to the nth focal zone, is constructed using a single set of delays (one delay per element). The set is 
obtained using Q, which is calculated for a single, representative point within the focal zone. 

In this section, we have formulated the manner in which the beamformer combines the signals received in the 
transducer elements, { ( Pm('t)} m —-M into a beamformed signal, a), through the process of dynamic focusing. 
The dynamically focused, beamformed signal, may now be used in order to generate a single image line. The 
beamformer performs the computation formulated in (|6j) in the digital domain, using samples obtained from 
each of the transducer elements at the Nyquist rate. Our goal is to retrieve a set of parameters from which the 
beamformed signal $(£; a) formulated in ([6| may be reconstructed, by sampling the received signals far below 
the Nyquist rate. We achieve this by exploiting an approximately FRI structure characterizing within 
the Xampling methods of Refs. [10] and [Tl] which are outlined in the next section. 

Throughout the rest of this paper, we will limit ourselves to the case of linear scan, in which all beams are 
parallel to the z axis. This is achieved by setting a = in The parameter S m now represents the distance 
between the mth receiver and the beam processed at the current cycle. Equation ^ then becomes: 

$(t;a = 0)= E y m Q(*+' v / * 2 + 4 (Wc) 2 )). (7) 



3. XAMPLING THE SIGNAL OBTAINED IN A SINGLE TRANSDUCER ELEMENT 

Up until this point, we have outlined the process in which a B-mode ultrasound image line is generated. The 
cycle begins by transmitting a modulated pulse along a narrow beam. The device then captures the intensity of 
echoes reflected along the beam by applying dynamically focused beamforming. 

Regarding the ultrasound signal detected by a single transducer element indexed m (denoted by tp m it)), Tur, 
Eldar and Friedman^ assume that it consists of a set of known-shape pulses, which result from reflections of the 
transmitted pulse by strong, macroscopic perturbations in the tissue. The signal can hence be approximated as 
an FRI of the form: 

L 

<Pm(t) = X! a l,mh(t - ti >m ), (8) 
1=1 

where hit) is a known-shape pulse, and there exists some r > such that h(t — t/ lTO ) = 0, Vt ^ [0,r), I = 1...L. 

The extremely short support of hit) implies that ip m (t) is of very wide band (the support of hit) is typically 
600 Ilsec , with r being typically 200 psec ). Classic Nyquist-Shannon sampling theorem thus forces standard ul- 
trasound devices to sample tp m (t) at a high rate (typically 20m Hz)- Nevertheless, we may easily observe that 
<Pm(t) actually has only 2L degrees of freedom (L being the number of macroscopic scatterers along the path of 



the transmitted pulse). The schemes introduced by Refs. [10| and 11 manage to exploit this property, enabling 



the reconstruction of ip m (t) from a much smaller number of samples (at least 2L per time interval r) 



In this section we outline the system suggested by Ref. 11 for Xampling the signal ip m (t)- This approach is the 



basis for our scheme, introduced in Section [4] We emphasize that both Refs. [TT] and [10] treat Xampling of a 
one dimensional signal which is received in a single transducer element. Moreover, they do not treat integration 
of samples, obtained from multiple elements, into a two-dimensional ultrasound image. Our novelty concerns a 
method for obtaining parameters of the beamformed signal, $(<; a = 0) (introduced in the previous section), 
which is directly related to the image line. We obtain these parameters from low rate samples of the individual 
signals tp m (t). By extracting the parameters of the beamformed signal we also cope with the noisy components, 
which induce erroneous results when attempting to reconstruct the individual signal (p m (t) , from its corresponding 
low rate samples. 

Let us denote by H(uj) the CTFT of the known-shape pulse h(t), and by </> m [fc] the kth Fourier coefficient of 
^ m (t)'s r-periodic extension. Further denote by k, a set of K consecutive indices for which H(u> = ^fk) ^ 0, 



Vfc £ k. Ref. 10 shows that, as long as K > 2L and the unknown time delays are distinct, i.e. U m ^ tj. m , 
Vi ^ j, one may accurately estimate <fi m (t), from the set {4> m [k]} k( z K - Gedalyahu, Tur and Eldai^ suggest a 
practical approach for obtaining the set {0m[fc]}fc £K , involving a bank of modulators and integrators. Referring 
to Figure [2] after having chosen the set of indices k, we construct p> K branches, each comprising a modulating 
kernel and an integrator. We then set the modulating kernels to be: 

M*)=E fce «(^e- J '^ fet ), 9 = 1,2,..., p. (9) 



Ref. 11 proves that the following relation holds: 

c = S(j), (10) 

where S is a p x K matrix with s q _k as its (q, k) element, c denotes the length-p sample vector with the output 
of the qth branch as its qth element, and <p denotes the length- K vector with the Fourier coefficient 4> m [k] as its 
fcth element. As long as S has full column rank, we can recover <\> from the samples by 4> = S^c. 

Denote by H the K x K diagonal matrix with fcth entry H(u> — —k), fee K, and by V (t) the K x L matrix 
with (k,l) element e _j22 ? Lfe *'' m , where t = {t± lTn , tL, m } is the vector of unknown pulse delays received in the 
individual transducer element. In addition denote by a the length-L vector whose Zth element is a^ m . Then: 



<i) = HV (t) a. 



(11) 
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Figure 2: Gedalyahu, Tur and Eldar's^ multichannel Xampling scheme. The resulting samples {c q } p =1 
constitute a mixture of Fourier coefficients corresponding to tp m (t)'s r-periodic extension. The coefficients' 
indices belong to the set k. Figure reprinted from Ref. 
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The matrix H is invertible by construction. Generating the length- if vector y by left multiplying by H 1 , we 
have: 

y = V(t)a, (12) 

which is a standard problem of finding frequencies and amplitudes of a sum of L cisoids (complex sinusoids). The 
time-delays {iz,m} ;=1 may be estimated using nonlinear techniques (e.g. annihilating filter,^ or matrix penciP^ 
methods). Having obtained the time delays, estimating the amplitudes {o,i,m\i—i is a linear problem, which may 
be easily solved using a least squares approach. 

In the next section we apply the scheme of Figure[2]on the beamformed signal $(t; a — 0), in order to reconstruct 
it from a small number of its samples. Recall that standard ultrasound devices digitally construct a = 0), 
after sampling the individual signals received in the transducer elements at the Nyquist rate. Since our goal is 
to break the Nyquist barrier, we bypass the actual construction of <I>(i; a = 0), by translating its Xampling to a 
scheme which may be applied directly on the analog signals {^m(0} m =-Ji/- 



4. GENERATING A 2D IMAGE BY XAMPLING THE BEAM FORMED SIGNAL 

At the basis of our approach is the assumption that the beamformed signal $(i;a = 0) maintains the FRI 
property which characterizes the signals from which it is constructed. This property was formulated in ([8]). 
More explicitly, <£>(£; a = 0) may be written in the following manner: 

L 

*(t;a = O)=5^6i/i(*-<0- ( 13 ) 

This claim requires justification, which will not be provided within the scope of this paper. The justification is 
concealed within the fact that the nonlinear scaling of ip m (t), formulated in ([7]), has little affect on the shape of 
the pulses which construct it, due to their extremely short support with respect to r. 

This notion is well demonstrated in Figure [3j A single pulse is transmitted along the narrow beam extending 
from the origin, along the z axis. Two scattering elements are illuminated by the beam (one is distanced l cm 
from the origin and the second is distanced 2 cm from the origin). As the pulse interacts with the elements, 
reflections are scattered and received in each of the 16 array elements. The traces well demonstrate the different 
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<I>(t;u=0) Generated by Summing the Distorted Signals 




Figure 3: Generation of = 0) as defined in ([7]), by summing distorted versions of the signals received 

by 16 elements. The setup is depicted in the top left: two scattering elements are positioned along the z axis, 
reflecting the pulse transmitted by an array of 16 elements. The 16 traces received by the elements are depicted 
in the top right. Notice that the reflected pulses are not aligned due to the different path traveled to each 
element. After applying the distortion suggested in Q, we obtain modified traces, in which the pulses appear 
aligned, at the cost of slight distortion to the shape of each replica (bottom left) . Finally, summing the distorted 
traces yields the signal <£> (t; a = 0) (bottom right). The latter may be approximated as a delayed sum of two 
known-shape pulses. 



pulse delays obtained in each of the elements. Simply summing the 16 traces (namely - beamforming with 
receive focus set to infinity) may yield certain SNR improvement, yet artifacts will be formed due to the fact 
that the corresponding pulses are not aligned. Using standard imaging techniques, these artifacts are manifested 
in a non-focused image. Yet they may have even more profound implications when attempting to Xample the 
beamformed signal. Let us now distort each of the received traces, parametrized by m, as defined in |7]): 

(M*) = <Pm Q (t + ^i 2 + 4 (WcfU . (14) 

Observing Figure [3j one may notice that the pulses are now aligned, although each may have undergone a slight 
distortion. Summing the distorted signals, we now obtain <i> (i; a — 0) which may evidently be regarded as FRI. 
If $ (t; a = 0) actually existed in the analog domain, then we could Xample it using the scheme described in 
Section [3j since it approximately satisfies (13). Namely, we could reconstruct it from a rather small subset of 
Fourier coefficients corresponding to its r-periodic extension. 

Let us now assume that $ (t; a = 0) existed in the analog domain, and feed it to the input of the scheme depicted 
in Figure [2j instead of the individual trace tp rn {t). Formulating the operation of a single branch, with $(i; a = 0) 



as its input, we have: 



M 
i=-M 



t+^t 2 +4{6 m /c)' 



dt. 



(15) 



After several algebraic manipulations, and exchange of variables inside the integral, ( 15 ) may be brought into a 
rather straightforward form: 



M fi r 

E \-J Sq,m{t)<Pm(t)dt 



1 



=-M 

where s qim (t) is defined as: 

Sq,m (^) 

for 1 < q < p, —M < m < M, and u (•) is the unit step function: 

, , f 1 i>0 
U( " )= else 



Ct ' 



E««.*«-^ fc(t -* (¥) '4«(*-i-i). 



The formulation of (16) may now be interpreted by the following algorithm, applied on the set {ip m (t)} 



M 



(16) 



(17) 



(18) 



-AT 



1. Having defined the full column-rank px K matrix S, generate the extended set of p x (2M + 1) modulating 
kernels s g>m (t), defined in (17). 

2. Modulate each of the analog signals received by the active elements (p m (t), m G {— M, ...,M}, using its 
corresponding, size p, set of kernels, {s 9jm (t)}^ =1 , yielding a corresponding size p set of coefficients: 



~ fo 8q,m(t)<Pm(t)dt q €{!,... ,p} 



(19) 



Applying the above step upon each of the signals </? m (t), m € {— M, ...,M} yields a p x (2M + 1) matrix 
of output samples. Note that the upper integration bound was modified to f, where: 



t = max 

mE{-M,...,M 



(20) 



3. Sum the samples along to, yielding the single, length-p vector, c, of which the qth element is c q , satisfying: 



M 



(21) 



m=-M 



4. Obtain the vector y by: 

y = H- 1 (St c ). (22) 

Since we began the derivation by injecting <!>(i; a — 0) into the scheme depicted in Figure [5J then the 
vector recovered by S^c, holds, as its fcth element, the fcth Fourier coefficient corresponding to the r- 
periodic extension of the beamformed signal $(t; a = 0). 



5. Solve the problem formulated in ( |12[ ), extracting the 2L unknowns of ( 13 ). Here V (t) is the K x L matrix 
with (k,l) element e -j2 ^ L ' £ * ! , t = {ii,...,^} is the vector of unknown pulse delays, and a is the length- 
L vector whose Zth element is 6;. The time delays are extracted by applying nonlinear techniques (e.g. 
annihilating filter, or matrix pencil methods). Estimating the amplitudes is then a linear problem. 



Summarizing the above, Figure [4] schematically depicts the suggested generalized Xampling scheme. The active 
array elements typically lie symmetrically with respect to the processed image line, resulting in symmetric kernels, 
with respect to m (i.e. s q - m (t) = Sg,m (£))■ Exploiting this symmetry enables reducing the total number of 
samples, by summing each pair ip m (t) and tp- m (i) prior to the modulation and integration branch. This concept 
is implemented in the figure. 

Recall that the transmitted (and reflected) pulses are modulated to a high frequency carrier (typically f c = 
^>MHz)i such that most of the pulse energy is concentrated far from the DC frequency. In order to avoid 
singularity of the matrix H it is necessary to pick the set k, such that the frequencies {ojk = ^^} fcgK are near 
27r/ c . This implies, that when implementing the Xampling scheme depicted in Figure [2] using real filters, the set 
K is doubled, so that if the index k is in the set k, we must also add — k to k. With the condition that p > |k| 
(so that S is full column rank), | • | denoting cardinality, we are forced to double the number of samples: for 
$(t = 0; a) having 2L degrees of freedom ({6;, ti} l=1 ), reconstruction requires a minimum of 2K samples within 



the interval [0,r), where K > 2L as before. Indeed, the actual ultrasound signal Xampled in Ref. 10 was first 
demodulated, yielding a "complex" analog signal. This practically means that the single, modulated real signal, 
was split into two base-band signals, each Xampled at least at the Rate of Innovation - thereby the number of 
samples obtained within the interval [0, r) satisfied p > 2K > AL. 

To conclude, we Xample the signals { ( fm(t)} m —-M m a manner which enables reconstruction of the beamformed 
signal $(t;a = 0), from which an image line may be directly generated: the individual delay U corresponds 
to a coordinate, located along the line normal to the set of active elements, at a distance ^ from its center. 
We may simply set the pixel associated with this coordinate an intensity proportional to 6;. Alternatively, we 
may generate a trace, by convolving the stream of pulses Yli=i bi5(t — t{) with the envelope of h(t), and then 
set the pixel intensities along the corresponding image line accordingly. The latter method was used for the 
simulations discussed in the next section, in order to obtain results comparable to those obtained by standard 
imaging techniques. 

For our simulations, we assumed L — 30 macroscopic scatterers along a single beam. Defining p > 1 to be 
the oversampling factor, we then need at least K = 2pL consecutive indices, for which H(uj = —k) ^ 0. We 
select the set of indices, k, such that — ^rk are near 27r/ c , where f c is the frequency of the carrier wave 
(approximately 5m Hz)- In order to obtain real kernels, we add the opposite indices to the set k, such that we 
finally have \k\ = ApL. We now construct a scheme with p channels per element, —M < m < M, such that 
p = |k| = ApL. We do this by choosing the matrix S in a very straightforward manner: 



±1 il 

2j 2j 



(23) 



where I is a | x | identity matrix. Referring to ([17]), our selection of S forms the following set of p x (2M + 1) 
kernels * 

m _J [i + (^) a ]«»(^(*-K¥) a ))«(*-l¥l) 

-[l + (^) 2 ]sin(^Vf (*-*(%0 2 ))«(*-lVO i + l<?<P ' { } 

which were utilized in the scheme depicted in Figure |4j 

We chose the matrix pencil method in order to estimate the set of delays from the estimated Fourier coefficients 
obtained using our scheme. One significant advantage of the matrix pencil scheme is that it provides the ability 
to estimate the actual number of delays concealed within the noisy data, based on a SVD decomposition process. 
Recall that in real imaging we have no prior knowledge regarding the number of reflecting elements aligned along 
a single image line. Nevertheless, when we design the matrix pencil we must determine the pencil parameter, 
which we shall denote by r\. The matrix pencil method assumes r\ to be greater than (or equal to) the number 
of complex exponentials comprising the estimated signal (in our context, this is equivalent to the number of 
delayed pulses, L). As a result, when designing the matrix pencil block we must first determine an upper bound 
on the number of reflected elements along a single image line and then set the pencil parameter r\ accordingly. 
This farther dictates a lower bound on the number of samples used for estimating the delays (denoted here by 
K). Summarizing the above constraints, the following must hold: 



L <rj < K - L, 



(25) 



where L is an estimated upper bound on the number of reflecting elements. The need to guarantee a non-empty 
interval [L,K — L], justifies the requirement that K > 2L. 




Figure 4: Generalized Xampling scheme yielding p coefficients, used for reconstructing Q(t; a = 0) from 2M + 1 
receiving elements. 



5. RESULTS 

In the following section, we examine the result of applying our suggested Xampling scheme upon actual raw RF 
ultrasound data. The data was acquired using a programmable imaging system (Model V-l-128, Verasonics, 
Inc., Redmond, WA), equipped with a 128-element 1-D linear transducer array (Model L7-4, Philips Healthcare, 
Bothell, WA). The imaging target was a commercial multi-purpose gray-scale phantom (Model 403GS LE, 
Gammex, Inc., Middleton, WI) including 0.1-mm nylon wires embedded in tissue mimicking material. We 
compare the following images: 

1. Standard Image - Ultrasound image generated using standard imaging technique: the frequency of the 
carrier wave was 5.142jvfH Z ; data was acquired at high rate (20mhz)', the number of active transducer 
elements used for focusing the beam at transmit (Tx) and at receive (Rx) was 16. Setting the maximum 
imaging depth to 7.88 cm , the number of samples used in order to generate a single image line was thus 
16 x 2048 = 32, 768. We used dynamic Rx focusing, comprising K = 100 focal zones. 

2. Xampled Image, Dynamic Rx Focusing - Ultrasound image generated using the exact scheme suggested in 
Section IU 

3. Xampled Image, Infinity Rx Focusing - Ultrasound image generated using the scheme suggested in Section|4j 
forcing S m = 0, Vm e {—M, ...,M}. Observing Q, one can easily see that this degenerate implementation 
forces Rx focusing to be infinity: $ (t; a) is simply the sum of all signals obtained from the active elements, 
without applying any delay. 



The resulting images are arranged in Figures [5]|9j according to the following table: 



Figure 


Image Type 


Estimated No. of 


No. 




Elements / Oversampling 
Factor 


5 


Standard 


Irrelevant 


6 


Xampled (Dyn. Focus and oo Focus) 


30/1 


7 


Xampled (Dyn. Focus and oo Focus) 


30/2 


8 


Xampled (Dyn. Focus and oo Focus) 


30/3 


9 


Xampled (Dyn. Focus and oo Focus) 


30/4 



Table 1: Properties of images displayed in Figures |5|9| 



Figure 5: Standard image comprising 113 Image Lines, generated with 100 focal zones. 




Figure 7: 



Xampled Images with L 



= 30 and p 



= 2. Dynamic Focusing (left) and Focus at Infinity (right) 




Figure 8: 



Xampled Images with L 



= 30 and p 



= 3. Dynamic Focusing (left) and Focus at Infinity (right) 



Figure 9: Xampled Images with L = 30 and p = 4. Dynamic Focusing (left) and Focus at Infinity (right). 



At this point of our work we avoided deriving quantitative measures to the estimation quality; this is because 
we had no actual documentation regarding the true positions and intensities of the reflecting elements. Using 
SNR measurements with the "Standard Image" as reference may be very deceiving, since the latter is obviously 
noisy by itself - noise which was often eliminated by the Xampling scheme, due to the fact that it extracts 
isolated delta excitations. The reader may receive qualitative impressions, observing the similarity between the 
Xampled images and the "Standard" image, and the obvious difference between a dynamically focused image 
and its equivalent infinity focused one. 

In Table [2] we provide an estimate of: a) the number of samples required for generating a single image line 
using our suggested Xampling scheme (referring to L and p used in Figures [6|9| and using standard imaging 
techniques; b) the estimated computational cost of generating the same image line for both methods: 



Image Type 


Estimated 


Oversampling 


K 


Sampling Rate 


Cost 




No. of 


Factor 




[No. of Samples/ 


[MegaOps./Line] 




Reflectors 






Element /Line] 




Xampled 


30 


I 


60 


120 


0.43 


30 


2 


120 


240 


2.81 


30 


3 


180 


360 


9.06 


30 


4 


240 


480 


21.05 


Standard 


Irrelevant 


2048 


0.06 



Table 2: Sampling rates and computational costs estimated for Xampled imaging (using various over-sampling 
factors) and for standard imaging. The computational cost of a Xampled image line roughly grows like O (|k| 3 ). 

The computational cost required for the standard imaging process (assuming 2048 samples per image line, and 
16 active receivers used for beamforming) comprises 2048 x 15 Add operations, and additionally the cost of the 
Hilbert transform which is utilized for envelope detection (typically implemented using two FFT operations). 
Table[3]which appears in Appendix A details the blocks used for estimating the Xampling scheme computational 
costs. 



6. CONCLUSIONS AND FUTURE WORK 



This work focused on generalizing the Xampling method suggested in Ref. [TT] to an array of multiple receiving 
elements. At the heart of our generalization is the observation, that the dynamic focusing and the filtering part 
of the Xampling can be combined into a set of modulating kernels and performed directly on the analog signals. 
This, in turn, can be sampled at a rate way below the Nyquist rate. 

Preliminary tests on actual ultrasound data yield results which are quite similar to an image obtained using 
standard techniques, while reducing the sampling rate by a factor of 5 to 15. Apparently, this is achieved at 
the cost of increased computational effort. By reducing the sampling rate, we hope to simplify the front end 
hardware (mainly in terms of size and power consumption) while maintaining image quality. 
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APPENDIX A. ESTIMATION OF XAMPLING SCHEME COMPUTATIONAL COST 



Operation 


Details 


Cost/O (•) 


Sum Outputs 


Sum (c ix( 2M+i),2) 


(L - 1) X (2M + 1) 


Matrix Pencil 
Obtain {ti}f =1 


Yjfxl = A^ X pCp X x 


K x P 


SVD (Q|kx|k) 




Y %Kx±K = U| KxL S' LxL V' 1 ix i K 


\K x Lx L + Lx Lx ±K 


Y2 Kx i K = U2 KxL T,' LxL V"- Lx i K 


\K x Lx L + Lx Lx 


pinv (Y2 Kx i k s ) 












Least Squares 
Obtain {bi}f_ 1 


LkxL * M-KxL 


KxL 


pinv (V KxL ) 


K 3 + L 3 


t>Lxl = A-LxKCKxI 


KxL 



Table 3: Blocks used for calculating Xampling scheme computational cost using the matrix pencil method. 



REFERENCES 

[1] Szabo, T. L., "Diagnostics ultrasound imaging: Inside out," in [Academic Press Series in Biomedical Engi- 
neering), Bronzino, J., ed., Ch. 7, 10, Elsevier Academic Press, 200 Wheeler Road, 6th Floor, Burlington, 
MA 01803, USA, first ed. (2004). 

[2] Jensen, J. A., "Ultrasound imaging and its modeling," Topics in Applied Physics 84, 135 - 165 (2002). 

[3] Jensen, J. A., "Linear description of ultrasound imaging systems." Notes for the International Summer 
School on Advanced Ultrasound Imaging, Technical University of Denmark (1999). 

[4] Shannon, C. E., "Communication in the presence of noise," Proc. IRE 37, 10 - 21 (1949). 

[5] Vetterli, M., Marziliano, P., and Blu, T., "Sampling signals with finite rate of innovation," IEEE Transac- 
tions on Signal Processing 50, No. 6, 1417 - 1428 (2002). 

[6] Mishali, M., Eldar, Y. O, Dounaevsky, O., and Shoshan, E., "Xampling: Analog to digital at sub-nyquist 
rates," IET Journal of Circuits, Devices and Systems 5, Issue 1, 8 - 20 (2011). 

[7] Mishali, M. and Eldar, Y. O, "Sub-nyquist sampling: Bridging theory and practice," submitted (April 2010). 

[8] Mishali, M. and Eldar, Y. O, "From theory to practice: Sub-nyquist sampling of sparse wideband analog 
signals," IEEE Journal of Selected Topics on Signal Processing 4, No. 2, 375 - 391 (2010). 

[9] Matusiak, E. and Eldar, Y. O, "Sub-nyquist sampling of short pulses: Part i," submitted to IEEE Trans. 
Information Theory [Online] arXiv 1010.3132 (2010). 



[10] Tur, R., Eldar, Y. C, and Friedman, Z., "Innovation rate sampling of pulse streams with application to 

ultrasound imaging," to appear in IEEE Transactions on Signal Processing (2010). 
[11] Gedalyahu, K., Tur, R., and Eldar, Y. C, "Multichannel sampling of pulse streams at the rate of innovation," 

accepted to IEEE Trans, on Signal Processing (2010). 
[12] Blu, T., Dragotti, P. L., Vetterli, M., Marziliano, P., and Coulot, L., "Sparse sampling of signal innovations," 

IEEE Signal Process. Mag. 25, No. 2, 31 - 40 (2008). 
[13] Kinsler, L. E., Frey, A. R., Coppens, A. B., and Sanders, J. V., [Fundamentals of Acoustics), John Wiley 

and Sons, New York, third cd. (1982). 
[14] Stoica, P. and Moses, R., [Introduction to Spectral Analysis], Prentice-Hall, Englewood Cliffs, NJ (2000). 
[15] Sarkar, T. K. and Pereira, O., "Using the matrix pencil method to estimate the parameters of a sum of 

complex exponentials," IEEE Antennas and Propagation Magazine 37, No. 1. 



